#include "stdafx.h"
#include "algo.h"

double gass(a,d,x)
double a,d,x;
{
	double y;
	if (d<=0.0) d=1.0e-10;
	y=0.5+0.5*errf((x-a)/(sqrt(2.0)*d));
	return(y);
}

static double errf(x)
double x;
{
	double y;
	if (x>=0.0)
		y=gam2(0.5,x*x);
	else
		y=-gam2(0.5,x*x);
	return(y);
}

static double gam2(a,x)
double a,x;
{
	int n;
	double p,q,d,s,s1,p0,q0,p1,q1,qq;
	if ((a<=0.0)||(x<0.0))
	{
		if (a<=0.0) printf("err**a<=0!\n");
		if (x<0.0) printf("err**x<0!\n");
		return(-1.0);
	}
	if (x+1.0==1.0) return(0.0);
	if (x>1.0e+35) return(1.0);
	q=log(x);
	q=a*q;
	qq=exp(q);
	if (x<1.0+a)
	{
		p=a;
		d=1.0/a;
		s=d;
		for (n=1; n<=100; n++)
		{
			p=1.0+p;
			d=d*x/p;
			s=s+d;
			if (fabs(d)<fabs(s)*1.0e-07)
			{
				s=s*exp(-x)*qq/gam1(a);
				return(s);
			}
		}
	}
	else
	{
		s=1.0/x;
		p0=0.0;
		p1=1.0;
		q0=1.0;
		q1=x;
		for (n=1; n<=100; n++)
		{
			p0=p1+(n-a)*p0;
			q0=q1+(n-a)*q0;
			p=x*p0+n*p1;
			q=x*q0+n*q1;
			if (fabs(q)+1.0!=1.0)
			{
				s1=p/q;
				p1=p;
				q1=q;
				if (fabs((s1-s)/s1)<1.0e-07)
				{
					s=s1*exp(-x)*qq/gam1(a);
					return(1.0-s);
				}
				s=s1;
			}
			p1=p;
			q1=q;
		}
	}
	printf("a too large !\n");
	s=1.0-s*exp(-x)*qq/gam1(a);
	return(s);
}

static double gam1(x)
double x;
{
	int i;
	double y,t,s,u;
	static double a[11]= { 0.0000677106,-0.0003442342,
		0.0015397681,-0.0024467480,0.0109736958,
		-0.0002109075,0.0742379071,0.0815782188,
		0.4118402518,0.4227843370,1.0
	};
	if (x<=0.0)
	{
		printf("err**x<=0!\n");
		return(-1.0);
	}
	y=x;
	if (y<=1.0)
	{
		t=1.0/(y*(y+1.0));
		y=y+2.0;
	}
	else if (y<=2.0)
	{
		t=1.0/y;
		y=y+1.0;
	}
	else if (y<=3.0) t=1.0;
	else
	{
		t=1.0;
		while (y>3.0)
		{
			y=y-1.0;
			t=t*y;
		}
	}
	s=a[0];
	u=y-2.0;
	for (i=1; i<=10; i++)
		s=s*u+a[i];
	s=s*t;
	return(s);
}